Crystallization in the vicinity of dynamical arrest 
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Dating from experiments more than 20 years ago, it has been reahzed that the crystaUization 
of hard colloidal particles in the vicinity of dynamical arrest has several anomalies, that render 
the conventional nucleation and growth model inappropriate. Subsequently, key researchers have 
shown the influence of gravity. Here we show that a simple lattice model can capture most of the 
phenomena associated with such systems. In particular, the model reproduces not only characteristic 
signatures of glass- forming systems, but also the interplay between quasi arrested dynamics and 
crystal nucleation. 



I. INTRODUCTION 

Dynamical arrest, that process in which many parti- 
cles dramatically slow in a concerted manner, has been 
the focus of considerable attention for some years [1-4] . 
Packing-induced arrest, usually associated with repul- 
sive interactions between the particles, occurs because, 
as density increases, the amount of space available to a 
typical particle becomes so small that it becomes effec- 
tively trapped by its neighbors [5, 6]. Even near arrest, 
though, rare particles in the system may, as a result of 
fluctuations, have somewhat more empty space immedi- 
ately available to them. If these rare spaces can be used 
for motion, and the empty space generated by the mo- 
tion passes throughout the system, then motion on long 
length scales can be generated. Thus, there have been 
efforts to describe dynamics using dynamically available 
volume, starting many years ago [7], and more recently a 
growing understanding of the connection between avail- 
able space and dynamics has begun to emerge [3, 8]. 

However, there are problems related to, but comprising 
more than the issue of arrest itself, whose understanding 
would have broader scientific impact across many disci- 
plines. For example, an understanding of the manner 
in which ordered structures grow in the vicinity of dy- 
namical arrest is one of the more pressing and impor- 
tant questions in the modern condensed matter science. 
In the early stages, nano-science greatly focusscd on the 
creation of increasingly small and more functional parti- 
cles, and devoted proportionately less attention to ratio- 
nal approaches to assembling structures from them. Now 
it is increasingly realized that useful devices will require 
us to fabricate ordered structures from these particles, 
and in the longer term it will be necessary to approach 
these questions in a fundamental manner. Research in 
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arenas from photonic crystals [9] to protein crystalliza- 
tion [2, 10, 11] has been hampered by the fact that high 
quality ordered structures are hard to make, and may 
require specific, complex, and expensive methodologies. 
For the vast majority of systems of practical interest, or- 
dered structures compete with dynamical arrest leading 
to micro-crystallites in a matrix of glass, or partially or- 
dered materials with poor coherence, involving defects, 
dislocations, missing layers. Usually these phenomena 
degrade the functional properties of the crystal, where 
they relate to technological issues [12]. 

Hard particles should be somewhat simpler to under- 
stand. However, experimental studies of ordering kinetics 
of such particles, though began some years ago [13, 14], 
are less understood than one might expect [15], given the 
importance of the questions. The deficit of data is be- 
ing rectified [16-20], but understanding of the existing 
results has been slow to develop, and many key issues 
are still not agreed within the community [14, 15, 20]. 

The experimental results, mostly based on the model 
system PMMA and cis-dccalin, are believed to be generic 
and are, in summary form, as follows. As a function of 
increasing particle density, the hard-particle system ex- 
hibits first a single-phase fluid and then fluid-crystal co- 
existence. Above the volume fraction ~ 0.545 the pure 
crystal is the equilibrium phase. This crystal phase grows 
via homogeneous nucleation up to (pg ~ 0.58, but beyond 
this there appears to be a sharp cross-over to (very slow) 
heterogeneous nucleation, mainly from the edge of the 
vessel, and free surface [13]. Contrary to normal expec- 
tation, the size of the crystallites formed decreases as 
one goes deeper into the crystal region (higher volume 
fraction) and one approaches arrest. This sharp phe- 
nomenon has been understood by many authors to be 
itself dynamical arrest, and theoretical treatments of it 
have been generally considered successful [5, 21] though 
these observations are complicated by new findings. In- 
deed, gravity also seems to play a role in hindering crys- 
tallization [14, 17]. Moreover, recent molecular dynam- 
ics simulations show that even in the absence of diffusion 
crystallization is possible in bulk homogeneous monodis- 
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perse hard-spheres [22]. 

In our approach we would hke to probe directly the 
emergence of dynamically arrested substances in a sys- 
tem that is also able to crystallize. Our aim, therefore, 
is to deal with both arrest and crystallization in a fully 
consistent manner and establish the connection between 
the two. We would also like to see how ordering phe- 
nomena are controlled by the numbers and heights of 
barriers to particle movement (on all length scales) that 
originate from the caging and to understand how (rare) 
empty space is managed in the process of forming ordered 
structures. 



II. THE MODEL 

Our model has been first introduced in Ref. [23] to in- 
vestigate the dynamics of glass-forming systems in the 
presence of crystallization. The model represents two 
phenomena typical of extended particles in a crowded en- 
vironment, and we shall explain it with the help of Fig. 1. 
The first aspect is that particles will tend to move from 
more to less crowded portions of space (see Fig. 1(a)) in 
order to mimimize their local free energy. The second 
aspect is that they may have to wait for sufficient kinetic 
energy to cross a local barrier, or wait for some (micro- 
scopic) time period for a neighboring particle to move 
aside in order to access a nearby space. Both of these 
effects can be represented by a barrier in the local free 
energy (Fig. 1(b)). The aim is then to understand how 
these microscopic energy scales lead to long length and 
time scale behaviors. 

First, we define a Hamiltonian which makes the 
more crowded environments energetically relatively un- 
favourable: 

V 

H^VrY, - cr) e{n, ~ cr). (1) 

Here cr is the maximum number of nearest neighbors 
that may surround a particle without it incurring a cost, 
Uj the number of nearest neighbors of the particle at the 
j-th site and Vr is the strength of the repulsive inter- 
action, 9{x) is the Heaviside function and V the total 
number of sites (volume) . This Hamiltonian can be seen 
as an extension of the Biroli-Mezard model, where point 
particles on a lattice are forced to keep some exclusion 
volume around them [24]. Such a type of interaction 
mimics, on a lattice, the behavior of hard spheres and 
therefore originates a crystal phase. In this paper we 
consider a soft repulsion to be able to mimic slightly soft 
spheres on a lattice [25] . 

The second phenomenon we represent is the so called 
local caging effect, in which the locations of the neigh- 
bors of a particle lead to a local free energy barrier to 
access a neighboring empty space. It has been shown by 
experiments [26] , and also reproduced by continuum sim- 
ulations [27], that dense colloids are characterized by this 



phenomenon. A particle spends most of the time rattling 
inside the cage formed by its neighbors and occasionally, 
if a fluctuation of the neighbors opens up the cage, it 
makes a longer movement and starts rattling again in 
another cage. As schematically shown in Fig. 1(a), the 
particle is caged by several neighbors in such a way that 
the illustrated movement is only possible if the original 
cage of particles fluctuates and opens sufficiently to pro- 
vide an exit path. The relative unlikelihood of this kind 
of fluctuations originates in the fact that there are rela- 
tively few such configurations. The local free energy will 
therefore reflect this in having a barrier between these 
two adjacent local minima (Fig. 1(b)). The barrier en- 
ergy scale is represented by a rate constant for single 
particle motion [4, 28, 29] which can be implemented in 
a Monte Carlo scheme. We therefore define the following 
kinetic rule. A particle can move from a site i to one of 
its empty nearest neighboring sites j only if: 

(a) the sum of its nearest and next-nearest neighbors 
is not larger than a fixed parameter ck] 

(b) the movement is reversible, i.e. the particle in j 
can go back to site i without breaking the rule (a) . 

This rule takes into consideration the intuition that the 
number of particles involved in the local cage is larger 
than the one which controls the amount of space avail- 
able for intra-cagc motion. Thus, we assume that the 
arrangement of particles involved in a local cage extends 
up to the next-nearest neighbors. Moreover, because of 
the discrete nature of the model, a wider range of values 
of ck allows a more satisfactory fine tuning of the kinetic 
constraint. 

This kinetic rule has been inspired by the Kob- 
Anderscn model [28], which reproduces many signatures 
of systems close to dynamical arrest, including blocked 
non-ergodic states and dynamical heterogeneities [3, 30]. 
Moreover, it has been understood that empty spaces of 
a highly dense lattice model can be categorized accord- 
ing to their role in the dynamics. Holes are defined as 
empty sites which a particle can move into and constitute 
the dynamically available volume of the system [3]. In a 
dense configuration, most holes will be only involved in 
local movements. However, for some of the holes there 
may be at least one sequence of movements which allow 
the hole to move almost all the particles in the system. 
Such holes are named connected holes and play a funda- 
mental role in determining the diffusivity of the system 
[3, 31]. An important characteristic of this class of ki- 
netically constrained models is the presence, in addition 
to the local cages, of large cages, i.e. extended closed 
arragements of particles which prohibit every particle in- 
side the cage to move outside. Only from rearrangements 
of the particles on the outside those cages can be broken 
[8]. The presence of extended cages in our model has 
certainly an effect on the process of crystallization that 
we want to investigate. 

The two previous features of the model can be imple- 
mented in a Monte Carlo scheme by defining the total 
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transition rate probability Pi^j of a particle going from 
site i to site j as the product of a kinetic term and an 
energy term: 

The energy term is the usual Metropolis rule 

P,^, = min{e-^^, 1}, (3) 

where AE is the difference in energy of the two states ac- 
cording to the Hamiltonian (1). The kinetic term, which 
implements the constraint we have discussed above, is 

Ki^j = 9{cK - ni)B{cK - rij), (4) 

where rii and Uj are the sums of nearest and next-nearest 
neighbors of the site i and j, respectively, and 9{x) is the 
Heaviside function, with the convention 6(0) = 0. The 
quantity T = f3~^ thereby represents the effective height 
of the local barrier. 




(a) (b) 



FIG. 1. (a) Schematic representation of the caging phe- 
nomenon. To go from a crowded (a; a) to a less crowded 
position (xb), the particle has to overcome the barrier (b) 
due to particles in the immediate vicinity (such as C and D) . 
Cage escape rates, that generally depend on the type of inter- 
action and particle density, are represented by kinetic rates in 
Kob- Andersen models [28]. AE is the difference in local free 
energy between cages. 



III. PROPERTIES OF THE MODEL 

In this paper, we want to focus on aspects related to 
crystallization and nucleation close to arrest. Evidently 
such a model has an extensive parameter space, which 
is indeed likely to represent a wide range of physically 
realizable situations. However, here we wish to explore a 
scenario that is close to the widely studied hard sphere 
case. Choosing cr = 3, Vr = 1, it is possible to cal- 
culate an equilibrium phase diagram which presents the 
required characteristics. A description of the equilibrium 
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FIG. 2. Section of the equilibrium phase diagram of the model 
for Cfl = 3 on a cubic lattice [25] , F=fluid, C=crystal. In the 
crystal phase for T < 0.4 and p > 2/3, signatures of a glass 
state (g) are observed for ck ~ 10 [23]. 



phase diagram [25] for such parameter values, as well 
as a detailed study of the dynamical properties [23] of 
this model have been given elsewhere. For the reader's 
convenience, we briefly summarize those results. Fig. 2 
reproduces an interesting section of the phase diagram. 
In that section, when the temperature is low enough, the 
phase diagram displays the sequence of phases (for in- 
creasing density: liquid, liquid-crystal coexistence, crys- 
tal), driven by repulsive short-ranged interaction, as in 
colloids. For cr = 3, the crystal of our model consists 
of double diagonal layers alternated with single empty 
diagonal layers and therefore it has a periodicity of \/3 
lattice steps. 

The origin of the fluid-crystal phase transition in this 
model can be explained by an entropic argument. The 
short-ranged repulsion represented by the Hamiltonian 
(1) has the effect of imposing some excluded volume in 
the system. Therefore, in the supercooled region the 
ordered crystalline structure is entropically favoured re- 
spect to the disordered fluid phase, due to the more ef- 
ficient use of the free spaces, as in hard-sphere systems 
[32]. The onset of crystallization in the supercooled re- 
gion is associated with the appearance of local assemblies 
of particles with the correct crystal structure. Such as- 
semblies are thermodynamically stable and we will refer 
to them as crystal nuclei. 

Inside the crystal phase, when we set ck ~ 10, an ap- 
parent dynamical arrest for density p > 0.66 and T < 0.4 
is observed (schematically indicated in Fig. 2). In this re- 
gion, several signatures can be found, such as extremely 
slow energy relaxation and non-exponential slowing phe- 
nomena, with good fitting of the Kohlrausch- Williams- 
Watts law. The calculation of the Kauzmann tempera- 
ture yields Tk = 0.42 [23]. 

This glassy behavior is determined by the presence of 
the kinetic rule, as in its absence the system quickly crys- 
tallizes [25]. The role of the parameter ck, in fact, is 
mainly to control the barriers and consequently both dy- 
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FIG. 3. Scheme of sequential growth of the crystal for cu = 3. 
In order to insert the particle marked in green, the kinetic 
parameter ck has to be larger or equal to 6. 
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FIG. 4. Energy evolution in number of lattice sweeps of a 
simulation a,t ck = 6, T — 0.4 and p — 0.64. This energy be- 
havior seems to indicate a steady state. Actually, the config- 
urations are disordered, but the equilibrium state is a crystal 
here. Therefore, the system is in a long-lived non-equilibrium 
state. 



namics and kinetics. Firstly, the value of ck is greatly 
involved in the detailed means by which crystals are 
formed. The sequential nature of the growing of layers 
or the ease with which one can fill a defect are features 
mainly governed by ck- Let us consider the crystal for 
the case cr = 3 and Vr = 1 at T = (i.e. in the hard 
sphere case) . It is interesting to note that a small value of 
Ck can determine quite strictly the mechanism by which 
the crystal can form. The minimum number of neigh- 
bors (i.e. the sum of nearest- and next-ncarest-nciglibors) 
which is necessary to fill a point defect into an isolated 
single diagonal layer of particles is only 3, but in order 
to fill a defect in a double diagonal layer, the number is 
9. A realistic scheme for growing a crystal at high den- 
sities is the case of a sequential growth on the border of 
an incomplete layer which is adjacent to a full layer (see 
Fig. 3). In this case, the lowest value of ck ior growing 
the crystal is 6. However, ck = G still appears to be a too 
low value to show an appreciable degree of diffusivity. As 
an example, Fig. 4 shows the energy evolution of a long 
simulation for p = 0.64 at T = 0.4. According to Fig. 2, 
this point is close to the left border of the crystal phase, 
but no crystalline state is observed. On the contrary, it 
is worth to remark again that we want to reproduce, as 
well as it is possible in a lattice model, the experiments 
which display classical crystallization and homogeneous 
nucleation as well as slow dynamics and arrest [13]. In 
that respect, the behavior at c/f = 10 is more interesting. 
Here two distinct dynamical regimes (crystallization and 
arrest) are easily identifiable. This observation suggests 
that in the process of nucleation the sequential growth 
is probably not enough to guarantee the crystallization. 
The possibility of filling point defects and, more gener- 
ally, of growing layers without a prescribed order plays 
a crucial role. Therefore, in the following results we are 
going to study extensively the model defined by ca' = 10. 



It is of some interest to understand the onset of this 
apparent dynamical arrest transition as temperature is 
lowered, and density increased, and the simplicity of the 
model makes it possible to dissect the mechanisms. For 
such models, changes in the nature of the dynamical pro- 
cesses (essentially the breakdown of the Stokes-Einstein 
relation) have been understood in terms of a dramatic 
reduction of the number of easy pathways in the phase 
space [3]. In turn this has been related to the manner in 
which excess space is organized, connected empty space 
allowing for long range dynamical process, and diffusion. 
Thus, if one creates an isoenergetic phase space (Kob- 
Andersen model) from the present model, by removing 
the energy term, the resulting model of the homogeneous 
fiuid exhibits a violation of the Stokes-Einstein relation. 
When the energy terms are restored, the number of free 
motions in the system is reduced yet further and, just 
beyond the Stokes-Einstein violation, the system appears 
to arrest. Nevertheless, one should not assume that this 
is necessarily a true glass transition, for the divergence 
of characteristic time is obtained by extrapolation, and 
there remains slow, but visible motions well beyond the 
apparent arrest. However, we can say with certainty that 
in this region the system becomes sub-diffusive, and the 
slowing is so dramatic that it is hardly possible to discern 
the difference from a real glass [23]. 

We observe that, for sufficiently high barrier heights, 
this phase diagram is similar to that shown in many 
experimental studies of the dense hard sphere system 
[13, 33]. Experiments also show that the kinetics of 
crystallization becomes different in approaching the glass 
transition. In the next section we seek to simulate those 
experiments using our model, as a possibility to check 
the quality of the mathematical representation and to 
give new insight into the phenomena involved. 
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IV. KINETICS OF CRYSTALLIZATION 

We are now in position to investigate in the model the 
phenomena associated with crystaUization. Here the aim 
is to understand to what degree the model reproduces 
nucleation and growth, and if it is possible to recognise 
dynamical arrest as in the experiments of colloidal parti- 
cles. In our analysis we will often refer to the two impor- 
tant characteristics of the model, namely the fluid-crystal 
phase transition generated by the short-ranged repulsion 
and the properties of the kinetic constraint, including the 
presence of extended cages in the system. Finally, we will 
try to understand if the model can even predict laws that 
could be checked experimentally. 



A. Configurations 

To begin with, we phenomenologically illustrate the 
kinetics of relaxation by showing configurations of dif- 
ferent densities as a function of time. After a period 
of equilibration, we see, depending on the system den- 
sity, fluid-crystal, crystal and "arrested" regimes as the 
density of the system is increased. This is shown in the 
sequence of figures (Fig. 5) where snapshots have been 
taken after a "long" period of time: meaning that all of 
the fast processes have ceased, and only slow aging re- 
mains. This representation may be directly related to 
the samples shown in the original papers [13, 33]. For 
the lowest density [p = 0.58), we sec. as expected from 
the equilibrium phase diagram, phase separation between 
fluid and crystal. Similarly, the next few densities quite 
rapidly equilibrate to the pure crystal state. However, 
with increasing density the system forms smaller crys- 
tallites that jam each others subsequent progress. Only 
after much larger times these crystallites coarsen to form 
a single crystal. Near the apparent arrest (p = 0.655) the 
size of the crystallites remains quite small. However, it 
must be pointed out that the samples at p = 0.655 and 
p ~ 0.665 will eventually crystallize after a longer time. 
As it was observed in the previous section, there are sev- 
eral indications that the hypothetical arrest should occur 
at p = 2/3, which is the density of the perfect crystal. 
Beyond p > 0.67 the system is slower and slower and the 
resulting states are more and more disordered at a fixed 
time. 

In Figure 6, another type of analysis is presented. We 
consider a glass state at high density (p = 0.69) and re- 
move a layer. In order to reduce as much as possible the 
finite size effect, wc consider an elongated sample. In 
this way, the removal of a layer on the small side pcr- 
turbates the system less, because it reduces the density 
decrease, so that the characteristics of the sample are not 
profoundly changed. Due to the excess of available empty 
space, the system partially crystallizes in the proximity 
of the emptied layer. However, the formed crystallites 
appear to stop growing at about t ~ 6 • 10^. In fact, the 
number of moving particles decreases sensibly and in the 



centre of the sample we observe a block of matter which 
appears quite resilient to movement. 

B. Analysis of the static structure factor 

In order to discuss more precisely these qualitative re- 
marks, we study the time evolution of the static struc- 
ture factor S'(fc), which can be regarded as an indication 
of the crystallinity of the system. We calculate the static 
structure factor of the configurations generated by our 
simulations using the formula 

5(k) = ^(pkP-k), (5) 

where the Fourier components pk of the density are cal- 
culated by the Fast Fourier Transform algorithm. 

A few comments should be made about the meaning of 
the structure factor in a lattice model. On a lattice, all 
the coordinates and the distances between pairs of points 
are discrete. For example, this implies that the number 
of pairs at short distances r is necessarily smaller than 
at larger distances. This causes a poor sampling aver- 
age in calculating the spherically-symmetric integrated 
structure factor S{k) at large k. Moreover, the behavior 
of S{k) in a crystal is quite different in a lattice model 
from a continuum model or experimental plots. Here the 
static structure factor of a perfect crystal consists of a 
single spike at the value of k corresponding to the peri- 
odicity of the crystal, and the peak has no width. On 
the contrary, in a real system molecules oscillate around 
their equilibrium positions and therefore slight changes 
in the measured periodicity occur. On the other hand, 
the structure factor of an imperfect crystal on a lattice 
presents some secondary peaks close to the main crystal 
peak, but this does not mean that they represent an ad- 
ditional periodicity in the sample. They can rather be 
considered as discrete effects of the width of the main 
peak, due to defects, dislocations or broken layers. In 
a real system, there is a much larger number of parti- 
cles and it is therefore highly probable that all possible 
slight modifications of the correct periodicity arc present. 
Therefore, we conclude that in a lattice model the height 
of the main peak of the structure factor quantifies more 
meaningfully the degree of order in the system than an 
integration over a number of close peaks. 

This is also the reason why we do not implement the 
definition of crystallinity recently used to analyze experi- 
mental data [18, 19], where the crystallinity is defined as 
the integral of the static structure factor over the peak as- 
sociated with the relevant periodicity. Moreover, it must 
be stressed that here we are not interested in the struc- 
ture factor itself, but in the divergence of the crystallinity 
in approaching dynamical arrest. As in the vicinity of 
arrest the growth of crystal nuclei involves a very long 
sequence of correlated movements (due to the increas- 
ing density), the region of space involved in the crystal 
growth becomes larger and larger, so that the divergence 



(a)p = 0.580 



(b)p = 0.640 



(c)p = 0.645 




(d)p = 0.650 (e)p = 0.655 (f)p = 0.665 



FIG. 5. Simulation of the experiment by [13]. Configurations in the fluid state at a given density p are instantaneously quenched 
at T = 0.4, where the equilibrium state is crystalline. Snapshots of the system after a relaxation process of 3 • 10^ Monte Carlo 
steps are presented. Blue spheres represent particles that have moved in the last 10^ time steps. Black spheres are particles 
which have not recently moved. 



of the characteristic times of the process cannot depend 
very much on the local details of the formed ordered zones 
in the sample. Proper crystallinity measures the amount 
of order in the sample, whereas the peak of the struc- 
ture factor simply represents an estimate of the size of 
the perfect crystal. We assume that these two quantities 
diverge with the same law in approaching dynamical ar- 
rest. Finally, we point out that the choice of using the 
height of the peak of the structure factor has been made 
in other experiments described in the literature [34] . 

In Fig. 7, the highest peak of the static structure factor 
S{kmax) is plotted against time for low densities, i.e. val- 
ues of p for which the crystalline state is accessible within 
our simulation times. Comparing the plot with experi- 
mental results [34], analogies and differences emerge. In 
experiments, two stages are usually recognised: an initial 
one, characterized by an intensity growth between and 
i**, and a second process with linear or sublinear behav- 
ior. Here (Fig. 7), it seems more a three-step behavior: a 
regime of slowly increasing S{kmax) is followed by a rapid 
growth and a plateau revealing that the state is fully crys- 
talline, as it can be verified by looking at the final states. 
As shown in the two fitting examples, the first process 



can be roughly described by sub-linear power laws, and 
the subsequent growth by a strong power law t" with 
a 4. This scenario, nonetheless, is compatible with 
classical nucleation: crystallinity grows up to a critical 
value and then jumps to the full crystal. As the melting 
point is pm ~ 0.637, the two lowest densities correspond 
to the two phase region, but the final value of S{kjnax) is 
roughly the same as the one of higher densities, because 
the fluid fraction is extremely small. 

As in the experiment by Harland and van Megen [35] , 
we define the induction time Tind as the end of the nucle- 
ation process, and the crossover time Tcross as the time 
when the equilibrium state has been reached; the growth 
time is At = Tcross — Tind (sec also Fig. 8). Quite interest- 
ingly, the value of S{kmax) at the induction time of Fig. 7 
appears to be roughly the same, with a slow increase 
closer to arrest. This looks slightly different from exper- 
iments, where the crystallinity seems to increase more 
with the volume fraction. The explanation is probably 
related to the fact that the crystallinity, based on the 
integration of the structure factor around the interesting 
peak, is a measure of the amount of order in the system 
and not a quantity proportional to the size of the largest 
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FIG. 6. Simulation of heterogeneus nucleation. A layer of particle is removed from the top of an "arrested" configuration at 
density p = 0.69 (the resulting density after the removal is p = 0.683). Snapshots of size 12 x 12 x 96 are shown during the 
relaxation process. Due to the extra free space, the system partially crystallizes. After t « 6 • 10^, the number of movable 
(blue) particles appear to decrease quite sharply. Blue and black spheres have the same meaning as in Fig. 5. We always use 
periodic boundary conditions. 



crystallite. Thus, if the density of nuclei increases with 
0, it is possible that the size of the critical nucleus stays 
roughly constant, whereas the number of nuclei increases, 
so that the sample exhibits a larger crystallinity, but the 
peak intensity remains constant. This could also give 



some insight into the phenomenon showed by the exper- 
iments after a time of the order of days [36], when the 
number of crystallites increases with the concentration, 
but their size gets smaller and smaller, up to dynamical 
arrest. That time is well beyond the crossover time and 
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FIG. 7. Highest peaks of the structure factor versus time for 
several densities lower than the hypothetical arrest transition 
at p — 2/3. Starting from a fluid state equilibrated at high 
T, we fix the temperature at T = 0.4 (so that the equilibrium 
state is the crystal). As an example, for p = 0.655 fits for two 
stages of crystallization have been calculated. Power laws 
with exponents a — 0.79 and a — 3.71 describe first stage 
of nucleation and rapid growth to the crystal, respectively. 
Increasing the density, the system crystallizes less and less 
easily: both induction and crossover times increase. 



cannot be directly inspected by our simulations, because 
of a finite size effect which obscures coarsening involv- 
ing length scales usually much larger than the simulation 
samples. 
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FIG. 8. Example of a simulation for p = 0.655, T = 0.4. The 
plot shows how the induction times nnd and the crossover 
times Tcross can be determined. t'„j represents our alterna- 
tive definition (see text). 

Approaching the apparent arrest transition, the inten- 
sity peak at Tind starts increasing more rapidly. The 
reason is that when the typical cage size becomes com- 
parable to the critical nucleus size, nuclei can grow larger 
than the critical value without jumping to crystal forma- 
tion and growth. In other words, the nuclei get caged by 
the kinetic constraint of the model and cannot grow any 
more. 

The situation beyond the apparent arrest transition 
(Fig. 9) is also relatively clear, as there is no evidence 



FIG. 9. S{kmax) as a function of time for two densities larger 
than 2/3. Deeply into the glass region, the crystal is not 
observed within the simulation time. A very weak increase 
of ordered assemblies is observed, but the time scale is enor- 
mously slower than for p < 2/3 (note the linear scale of the 
vertical axis). 



of critical nucleus formation on the accessible time scale. 
The question arises as to whether the arrested phase is 
simply a continuation of very slow nucleation beyond 
times accessible to simulations/experiments. It is also 
possible that the absence of gravity prevents any type of 
sedimentation and so true arrest is never reached. 

In the plots of Figures 7 and 9, the structure factor was 
averaged over a number of different simulations. This has 
been done to reproduce the global effect which should be 
observed on larger scales, as in the experiments. How- 
ever, this procedure actually hides the behavior of the 
single simulation runs, which present an interesting be- 
havior. Thus, in Fig. 10 and II, S{kmax) is plotted for 
several samples at densities p = 0.645 and p = 0.655. 
It is quite evident that different runs have different in- 
duction times. This scenario is unexpected, because it 
is not observed in typical first order transition models, 
such as the Ising model. As one can notice by comparing 
the two plots, this phenomenon becomes more impor- 
tant for higher densities, i.e. the spreading of the induc- 
tion/crossover times increases with p. 

According to the classical nucleation theory, nuclei 
growth is only hindered by surface tension. As soon as 
a nucleus reaches the critical size, it grows quickly over- 
coming the unfavourable surface contribution [37]. Here 
the picture is quite different. Figures 10 and 11, rep- 
resenting the time evolution of the highest peak of the 
structure factor at different densities, show that a new 
process takes place. Such process can be easily inter- 
preted by the language elaborated in the past years for 
kinetically constrained models [3, 31]. Those models are 
characterized by the presence of large scale cages and a 
diverging dynamical correlation length, which prevents 
diffusion at large density. In our model, this process in- 
teracts with nucleation and growth, because in a dense 
system the average cage size increases, becoming compa- 
rable to the critical nucleus size. At this point, it becomes 
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very probable for a crystallite to get trapped into a cage, 
and so it cannot grow as in the classical picture. In 
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FIG. 10. Plots of the time evolution of the main peak of the 
structure factor for many simulations: 50 samples for p — 
0.645. 
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FIG. 11. Plots of the time evolution of the main peak of the 
structure factor for many simulations: 50 samples for p — 
0.655. 

Fig. 10 we can single out three main behaviors, (a) A set 
of runs do show a behavior which is classical: the peak 
of the structure factor reaches a critical value and then 
grows rapidly to a large crystal. In this case, the critical 
nucleus is not caged and therefore it can grow without 
constraints, (b) In many other cases, S{kmax) shows a 
quite erratic behavior. The system does crystallize, but 
a large amount of time is spent in a regime where there 
are critical nuclei which do not grow any further. The 
result is a considerable increase of the induction time. 
This process can be interpreted as follows. Many critical 
nuclei develop in cages, and shrink or disappear before 
the cage is broken. Only when a critical nucleus forms 
in a region where it is not caged, the system can develop 
nucleus growth, (c) The few cases in which the system 
never reaches the crystalline state are caused by the fact 
that there is not any connected hole allowing the correct 
set of movements to form a full crystal [31]. In other 
words, the connected hole density is so low that a frac- 
tion of the samples does not contain any connected hole 
useful to move particles into the structure of the equi- 



librium crystal. Finally, it is worth underlining that this 
caging process becomes more and more important ap- 
proaching arrest. The growth is more and more delayed 
and S{kmax) evolution is characterized by non-monotonic 
behavior (see for example Fig. 14). 



C. Characteristic times 

We focus now on the distribution of induction and 
crossover times on approaching arrest. Each simulation 
presents a quite sharp transition to the crystal. There- 
fore, except for a few particular cases, the single sample 
induction time is almost coincident with the crossover 
time. Thus, we give an alternative definition of induc- 
tion time (Tj'„j^), as the instant after which the growth of 
S{kmax) has changed its speed and the sample clearly 
evolves to the crystal without large backward moves 
(Fig. 8). 

In Figures 12 and 13, the distributions of the induction 
and the crossover times are illustrated for the densities 
p = 0.645 and p = 0.655. Given the limited number of 
samples (150), it is hard to make definitive statements. 
However, it is possible to notice a few interesting phe- 
nomena. For example, it is evident that the width of the 
distributions increases with p. More tentatively posed, 
such distributions do not appear to be gaussian and seem 
to be skewed towards higher times. Moreover, it is inter- 
esting to note that approaching arrest this asymmetry 
seems to increase slightly towards higher times. 

As we have seen up to here, the kinetic processes that 
are still within the crystal phase but very close to dy- 
namical arrest (see for example Fig. 14) are quite pecu- 
liar, sufficiently so that we can no longer clearly identify 
the conventional nucleation and growth regimes. Instead, 
order grows in the system via a series of quite sudden 
changes, followed by periods where the system appears to 
be almost immobile. Furthermore, quite unlike conven- 
tional crystallization, the peak of the structure factor is 
no longer monotonic, and order can at first increase, then 
diminish for periods of time. This is a consequence of the 
wrong type of order having been formed, leading to a sort 
of "blind alley" , from which the system ultimately has to 
reverse, and try again to find a more successful pathway. 
In contrast to the classical picture, in the vicinity of ar- 
rest many nuclei grow up to few lattice steps and then 
jam each other instead of growing further. This creates 
a long-lived state, in which most particles are blocked or 
move locally, especially in proximity to the interfaces be- 
tween nuclei. This phenomenon does not appear to be a 
finite size effect, because simulations at larger sizes show 
similar crystallite sizes. To our knowledge, this is the first 
time such an intermittent mechanism has been observed 
in systems ordering near arrest; we suggest that these 
processes may be quite general, being simply related to a 
hard repulsive interaction and the mechanism of caging. 

The induction and crossover time of colloidal particles 
have been reported in a few papers [18, 35, 38]. Harland 
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FIG. 12. Distribution of the induction times calculated on the time evolution of the main peak of the structure factor for 150 
independent runs each. The densities chosen are p = 0.645 (a) and p — 0.655 (b), respectively. 
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FIG. 13. Distribution of the crossover times calculated on the time evolution of the main peak of the structure factor for 150 
independent runs each. The densities chosen are p — 0.645 (a) and p — 0.655 (b), respectively. 
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FIG. 14. Example of the time evolution of a sample very close 
to dynamical arrest (p — 0.659). Peaks of S{k) are plotted 
against time. It is noticeable that the curve is non-monotonic. 



and van Megen [35] studied the same type of particles 
described in the experiment by Pusey and van Megen 



[13], with an effective radius equals R = 201 ± Inm, the 
polydispersity was s w 5% and the fusion and melting 
volume fractions were (/i/ = 0.494, (pm = 0.545 ± 0.003, 
respectively. Schope et al. [18] studied similar particles 
with a radius R = 320nm, polydispersity s = 4.8%, and 
4>f ~ 0.505, 0m = 0.538; these two last values are based 
on the measurements of Kofke et al. [39]. For the data 
by Schope et al., dynamical arrest is estimated to occur 
at = 0.575. We now make a direct comparison between 
the experimental data and the results from our model. In 
order to do that, it is necessary to translate the plot of 
our characteristic times so that the melting point of the 
model, 0m = 0.637 coincides with the melting point for 
hard spheres 0m = 0.545. Then, we shrink our data on 
the p axis in order to make the hypothetical arrest tran- 
sition of the model 0g = 2/3 coincident with the hard 
sphere one: 0g = 0.58. Finally, the characteristic times 
have been multiplied by 0.1, so that they can fit the range 
of the experimental data. This is acceptable, because 
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the simulation times do not have an absolute meaning, 
but they are related to real time through a multiplica- 
tive factor. Fig. 15 compares our rescaled induction 





Hcirliiiid. 
Heymaiiii 
Scliope 4.8 
Scliope 5.3 
model 


1 

- - + — 

- K- - 

-o- 


1 

G 








, ' , ' 


' - ' ' ' P 




_ .6' " ' - ' ' 

.-r-.7.-.-.-:.t5 














* 1 "-t 



0.55 0.56 0.57 0.58 





FIG. 15. Comparison of plots of the induction time Tind versus 
density for several experiments and results from our model. 
The experimental data are taken from the papers: [18, 35, 
38] (s = 4.8%), respectively. H. J. Schope also provided us 
with data for polydispersity s = 5.3%. The values of Ti„d 
for the model have been multiplied by 0.1, so that they can 
fit the range of the experimental data. Here we focus on the 
asymptotic behavior close to the arrest transition, so what 
we compare is only the functional form of the divergence, and 
not the exact values of the characteristic times. 



the crossover times is illustrated. Unfortunately, in both 
cases the picture is not very clear, and even the different 
experimental sets are quite heterogeneous. The experi- 
mental procedure partially affects the data and polydis- 
persity also has a significant effect, as can be understood 
by a simple comparison between the two data sets from 
Schope et ai. In spite of these weaknesses, a similar trend 
is recognisable at high (f) in both induction and crossover 
times. In particular, a comparison between the asymp- 
totic behavior of the model and the data from Harland 
and van Megen (the ones closest to arrest), is intrigu- 
ing. There appears to be a quite similar functional form 
in the last points, a sign that could reveal an analogous 
divergence. 
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FIG. 16. Comparison of plots of the crossover times versus 
density for several experiments and results from our model, 
as in Fig. 15. Here we focus on the asymptotic behavior close 
to the arrest transition. 

times with sets from several experiments. Looking at the 
model results, it is evident that induction times Tmd grow 
sharply as we approach the arrest transition. The rapid 
increase of the induction times shows that it becomes in- 
creasingly difficult to form a critical nucleus that would 
liberate enough new free volume around it to make the 
growth process more rapid. Even when the critical nu- 
cleus is formed, the process of rearranging the system 
around it becomes more difficult. Therefore, it is appar- 
ent that the classical picture of nucleation and growth 
is no longer applicable. In Fig. 16, the comparison of 



FIG. 17. Growth times Ar vs for the model and experi- 
ments as in Fig. 15. 

The duration of the growth regime (Ar) also increases 
as dynamical arrest is approached. In Fig. 17, we present 
a comparison of the behavior of Ar vs p between our 
model and experimental data as in Figures 15 and 16. 
Again, our results have been processed in order to make 
(\)ra and 0g of the model coincident with the experimental 
values. Quite remarkably, the asymptotic behavior of the 
model and the data from Harland and van Megen appear 
to be similar on approaching arrest. No theory for either 
of these laws exist as yet, though it is natural that, as 
overall system mobility vanishes, they would also acquire 
a related divergence. 



V. CONCLUSIONS 

We note that our simple lattice model has been able to 
reproduce a wide range of phenomena from real glasses 
and energy landscape models, including onset of collec- 
tive behavior, divergence of the characteristic times, and 
many properties of crystal nucleation in the vicinity of 
dynamical arrest. The reason of the good performance 
of the model can be interpreted as follows. Models based 
on kinetic constraints alone do create a complex effective 
free energy landscape in which, at sufficiently high den- 
sity, many movements are prohibited by infinite or large 



12 



barriers. Nonetheless, many dynamical pathways (each 
mediated by a connected hole, in the language of the sim- 
ple models [3, 8, 31]) involving long ranged transport still 
remain at fixed (zero) energy. In such models, true dy- 
namical arrest only occurs when the lattice is fully filled. 
This dramatic reduction of dynamical pathways induced 
by kinetic constraints certainly leads to dynamical slow- 
ing, but not to true glass behavior. If we now allow 
different local energies within different cages one obtains 
a complex energy landscape. Then, rare pathways that 
were formerly barrier-less remain "easy" , but acquire a 
multitude of smaller energy barriers. The accumulation 
of such bumps against the backdrop of a vanishing num- 
ber of easy pathways ultimately leads to interesting sin- 
gular behavior for the characteristic time, that is consid- 
ered to be truly representative of the glass state. This 
is the root of glass behavior, and its physical origins arc 
quite clear in our model. Similarly, the presence of some 
of these pathways allow for the formation of the crystal. 

The only property that is not included in the model 
is the effect of polydispersity, which will be examined in 
a future project. However, a recent numerical work [22] 
shows that the effect of polydispersity on the dynam- 
ics is quite small and crystallization can also occur in 



a monodisperse system of hard spheres. Quite interest- 
ingly, from the analysis of our model it emerges that dy- 
namical slowing and anomalous crystallization appear to 
be essentially driven by the interplay between the kinetic 
barrier and the underlining crystal phase, as illustrated 
above. 

It must be considered intriguing that the ingredients 
of the model arc already known in the literature but that 
they have not hitherto been combined in this way. Purely 
repulsive Hamiltonian lattice models without kinetic bar- 
riers [24, 40] appear not to yield a KWW characteristic 
time law, as in experiments and continuous simulations. 
On the other hand, purely kinetic models do not have 
a crystal phase [28]. Here we have a simple model that 
reproduces, in spite of the weaknesses implied by the dis- 
cretization of the space, several main effects associated 
with colloidal systems. Thereby it opens up the possibil- 
ity of a more transparent dialogue between experiment, 
simulation and theory. 
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